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Abstract 



This paper presents a new approach to accommodating large, fluctuating renewable 
infeeds in power systems, based on the concept of reserve policies. The approach can be 
viewed as a generalization of existing reserve modes, and uses robust optimization with 
recourse to determine operating rules for power system entities such as generators and 

■ storage units. These rules, or policies, set in advance how these entities are to respond 

■ to errors in the prediction of loads and renewable infeeds once their values are discov- 

■ ered. We show how affine policies can be chosen such that the power network constraints, 

■ namely matching supply and demand, respecting transmission line ratings, and the op- 
erating limits of power system entities, are satisfied for all possible realizations of the 
prediction errors. The errors may arise from any number of spatially- and/or temporally- 
correlated sources. AfRne policies are compared with existing reserve operation modes 

0^ . under standard modelling assumptions, and operating cost reductions are reported for a 

multi-day benchmark study featuring a poorly-predicted wind infeed. Efficient prices for 
such "policy-based reserves" are derived, and we propose new reserve products that could 
be traded on electricity markets. 

(N 



1 Introduction 

A key challenge in incorporating highly variable intermittent renewable energy sources into 
power systems is the need to maintain system integrity while making best use of the energy they 
provide, which comes at zero marginal cost. It is widely agreed that in the next few decades, 
as the share of wind power becomes very large, current techniques for accommodating wind 
variability will become sufficiently expensive that alternatives will be sought [I112]- 

Running power systems with very high wind penetration and without excessive frequency 
control costs, or resorting to curtailment of renewable output, requires intelligent use of the 
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best available forecasts, at all times. In particular, any successful method for dealing with the 
high variability of renewables on intraday timescales (the only timescales over which prediction 
errors are reasonably small [H]) will require the following: 

1. A forecast of future intermittent energy injections available at the time when control 
decisions are made. 

2. Rules for acting on errors in this forecast when they are discovered. 

3. Forecast error probability distributions and their correlations in both time and space 
over the grid. 

Theoretical attention to these points has grown in the last few years as the share of wind 
power in several countries has grown [3HS]- The incorporation of uncertain wind inputs has 
typically been tackled using a scenario-based stochastic approach. For example, in Morales 
et al. [B] a unit commitment integer programming problem was solved first, and then in a 
second stage reserve margins were selected based on the requirement that the actual reserve 
deployment be feasible for all the scenarios considered. These ideas have since been developed 
to provide probabilistic guarantees on transmission constraint satisfaction using a limited 
number of scenarios [8]. 

In this paper, we also consider how reserves could be operated more efficiently on a daily 
timescale around any unit commitment decisions that have already been made. In contrast 
with existing literature, we use a robustness formulation of the problem and show that the 
use of affine policies can reduce expected operating costs in the presence of extreme energy 
fluctuations. A policy is in this case a rule, agreed in advance, governing how individual power 
system entities will respond to prediction errors affecting power system operation. It therefore 
constitutes a reserve mechanism, a concept that represents the main contribution of this work. 

This work was inspired by results on disturbance feedback policies from the control litera- 
ture. However these are predated by the concept of linear decision rules (LDRs) from opera- 
tions research, where current states, possibly together with past data or future predictions, are 
combined linearly in order to make an operational decision. Early examples appeared in the 
1950s, when Holt et al. used quadratic programming to schedule workforces in factories under 
uncertain future workloads [9]. Further applications included macroeconomic planning |10] . 
management [11) . and inventory control [K]. Typically, though, LDRs were unable to deal 
rigorously with operating constraints, and were not studied in much detail after the 1970s [13] . 
In the last decade, however, advances in constrained robust optimization have led to a revival 
of LDRs, for solving optimization problems where the minimizer is allowed to be a function 
of the data uncertainty |13H15j . This has led to some new applications, for example linear 
decision rules in portfolio optimization |16] . 



These solution methods were also shown to be applicable to robust predictive control 
[T7HT9] . where control policies with various dependences on the disturbance have been studied 
as a means of respecting state and input constraints under uncertain system dynamics. These 
control policies have since been applied in areas such as building control [20j . 

Although optimal short-term operation of power systems, including reserves, has been 
studied in various ways for decades as a variant of the standard optimal power flow problem 
|21j . recent advances in optimization over affine policies have not been exploited for real- 
time decision making in electricity provision under uncertainty. This application is interesting 
because policies may be an attractive way of using up-to-date forecasts in power systems. 
To this end, we present systematic ways of using future wind prediction error measurements 
to reduce the expected cost of reserve provision. This is made easier by the approximately 
linear transmission constraints present in power systems, namely matching consumption with 
production, and respecting line congestion limits. 

In Section [2] we show how generic state space models can be used to represent a wide range 
of power system entities. The electrical grid modelling assumptions are also given. In Section 
[3] the concept of policy-based reserves is described with reference to power system control. A 
finite-horizon optimization problem is defined and subsequently reformulated to facilitate a 
numerical solution, and an analysis of some properties of affine policies for power systems is 
given in Sectional In Section [5l we discuss optimal pricing and the practical implementation of 
policies. In Section[6]the potential benefits of policy-based reserves are reported in a numerical 
case study, and Section [7] concludes. Preliminary versions of some results appearing this paper 
have previously been reported in [22]. 

Notation: The use of (•)' denotes a vector or matrix transpose, [-Jfc denotes element A; of a 
vector, {X, Y) denotes the trace of product X'Y, and <^ denotes the Kronecker product. It 
denotes an identity matrix of size T. 

2 Power system model 

This paper considers the problem of optimal operation of an electrical network to satisfy loads 
in the presence of uncertainty. The uncertainty to be accommodated manifests itself in the 
form of random power infeed from renewables and fluctuating load requirements. We will 
choose operating rules that apply only for a finite time into the future, on the assumption 
that new rules will be determined before the present ones expire. This finite time, or planning 
horizon, is divided into T discrete time steps, corresponding to the Programme Time Units 
(PTUs), of length 5 minutes to 1 hour, over which electricity is traded on modern intra-day 
markets [23]. The length of time horizon considered relevant to this work is up to 24 hours, 
after which predictions of renewable infeed are assumed to become too poor to incorporate 



into sophisticated decision-making rules, and unit commitment decisions are not yet fixed. 
2.1 Participant model 

We consider the actions of Np generic entities, or participants, connected to a transmission 
grid. Each participant i, for example a generator, load, or storage unit, injects power into or 
extracts power from a fixed location on the network, in two forms (one or both of which may 
be present for a given participant): 

• Inelastic power flows, which cannot be influenced by control signals. 

• Elastic power flows, which are determined by the result of an optimization over possible 
control actions. 

Inelastic power flows 

The inelastic, or exogenous, injection or extraction of power for each participant i is modelled 
as Tj + Gi5, with positive values denoting a net power injection. Its two components are 
a nominal prediction G R-^ plus a linear function Gi G R-^^^*^ of entries of a random 
forecast error vector 6 G R^'^^, whose value is to be discovered in the future. It has the form 
5 = [6q, . . . , 6'j,_-^]', where each 6k G R^*. In other words, Ns is the number of elements in the 
disturbance vector at a given time, and this vector is mapped to the exogenous power flows 
in the system at that time. If the prediction error 5 is zero, then the net power injection at 
step j is simply [ri]j. 

The forecast error 6 is assumed to belong to a compact set A := {6 \ S5 < h} with /i G R"^, 
whose interior contains the origin. Although the error is random, the mean prediction error 
E[(5] and the covariance matrix E[(5(5'] G R^'''^^^^'^, which is symmetric, are assumed to be 
known. No other restrictions are placed on the probability distribution. 

Elastic power flows 

Elastic power flows are governed by a participant's dynamics in conjunction with some pre- 
defined costs. We describe them in state space as is common practice in systems and control 
( [24] . §2.1). At time k, each participant i has internal state G R"% where nj is the state 
dimension, and is governed by linear time-invariant dynamics, so that given an input at 
time k the state at time A; + 1 is given by a;^^;^ = Aixl. + Biul., where Ai G R"-»^"-» and 
Bi G R"'. The first element [x^i of the state vector is assumed to represent the current 
power injection at the relevant node of the transmission network, and other elements are used 
to model internal dynamics or memory of previous states. The scalar input u\ G R controls 
the net power injection of the participant at time k + 1. 



Assigning the current time the value A; = 0, a vector of future states for participant i, x* := 
[x\' . . . Xj^'Y S M"''^ can be written as a function of the input sequence u* := [uq' . . . Uj^_i]' G 
R-^ and the current state XqI 

x^ = Aix'o + Biu' , (1) 

where 
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The vector of outputs CjX* as seen by the network is just the power injected or consumed by 
the participant. Therefore each matrix Ci G M"^^"*^ selects only the first element of the state 
vector at each time, i.e. Ci = It ^ Ci, where Ci = [1 Oix(n,-i)]- 

Costs 

The function Jj : R"'-^ x MJ' — )• M is used to define costs for the states and inputs along the 
time horizon. It is assumed that Jj is a quadratic function 

J,(x^ u*) := /f x^ + ^x^'i/f X* + /r'u* + l^^'Hru^ + c, , (2) 

where the Hessian matrices fff and are assumed to be positive semi-definite, in order for 
the optimization problem defined in Section [3] to be convex. Linear components /f and are 
of the form Itxi ® fi and Itxi ^ ff respectively, where /f G M"* and /f € M. Similarly the 
quadratic components are given by := It ® Hf and Hf := It H^, where Hf G M"*^"-* 
and Hf" G M. Coupling costs between time steps can be represented by augmenting the state 
vector to include a memory of prior states in the state vector. Constant := Tq allows for 
a constant stage cost Cj. 

Constraints 

The set Zi consists of permissible combinations of state and input sequences x* and u* for 
participant i, which may in some cases be additionally constrained by 5. It is a compact set 
defined by k linear inequalities (i.e. a polytope) and takes the form 

Zi := < 

where Ti G R'»x"»^, Ui G M^'^^, G M'»x^«'^ and Wi G Of course, x' and u* are related 
by the dynamic equation ([T]), and Tj, Ui, Vi, and Wi may also depend on the current state Xq; 



X 



r,x* + Ku* + Vi5 < Wi 



(3) 



these dependences are left out of the notation above for clarity. Note that since x* and u* are 
trajectories rather than state or input vectors corresponding to a single time, a wide range 
of constraints coupling states and inputs can be modelled. For example, ramp rates may be 
imposed on generators, and empty/full constraints may be imposed on storage units. Binary 
variables, allowing for example start-up and shut-down constraints, are not considered in this 
paper. 

Usually Vi = 0, unless the uncertainty feeds directly into the participant's operating con- 
straints. An example of this would be a curtailable wind farm whose maximum power avail- 
ability at any given time is uncertain, and whose power output could be varied at any time 
between zero and this upper limit. In this paper, though, wind curtailment is assumed to be 
undesirable. 

2.2 Network model 

The network model is a standard linearized approximation of a high-voltage transmission 
grid [25J, in which lines are lossless, voltage magnitudes are constant, and line flows are 
proportional to the phase differences (assumed to be small) between nodal voltages. Let 
each participant be connected to one of network nodes, and let L be the number of lines 
connecting these nodes. 

The network imposes two constraints on power system operation. The first is that the net 
power injection, comprising the sum of inelastic flows [rj + Gi6]k and elastic flows [CjX*]^, has 
to be zero at all times fc = 1, . . . , T. This can be modelled using an equality constraint with 
T rows: 

Np 

+ G^6 + Ci^') = . (4) 

1=1 

The second constraint is that line currents cannot exceed the respective line ratings any- 
where on the network, at any time. Under the preceding assumptions, this constraint is linear 
in the power injections, as long as the net power injection into the network is zero (i.e. condition 
(jl]) holds) |25j . Hence an additional vector inequality constraint 

Np 

Y,Ti{ri + G^5 + a^')<p. (5) 

i=l 

This has 2LT rows, one for each flow direction, for each line, at each time. Each matrix 
Ti G M^^-^^^ maps the power output of the node to which participant i is attached to contri- 
butions to line flows. Each Fj can be constructed from the network line impedances using the 
derivation of equation (in.4) in [25], since we model line constraints as limits on the phase 
angle differences between adjacent transmission buses. 



3 Choosing optimal reserve policies 



Current reserve mechanisms use a cascaded loop structure, where the fastest (primary) con- 
troller stabilizes the grid frequency, a minutes-scale (secondary) controller corrects it back to 
its reference, and slower, separately-purchased tertiary reserves redispatch generators in order 
to free up the margins within which the faster control operates [23]. This reserve action is 
only a real-time response to the error as it unfolds, and makes no systematic use of what is 
expected to happen in future. Intraday electricity markets in many countries are currently 
experiencing dramatically increasing trade volumes [26], and this increase can be seen as an 
attempt to adjust the short-term economic operation of the power system in the light of new 
forecast information. Different countries operate these markets in diverse ways; an overview 
of the various mechanisms used to acquire reserves and short-term power commitments in 
European countries can be found in [27J. Such trading actions take little systematic account 
of time-coupled costs and constraints imposed on the market participants. 

In this section we describe a more systematic predictive mechanism that explicitly takes 
account of short-term uncertainties with the aim of reducing the expected running costs of 
the power system over the time horizon. 

3.1 Finite horizon optimization 

Consider the problem of minimizing expected running costs Ylf=i ^[Jii^^^ u*)] over a horizon 
of length T, subject to the local constraints ^ and network constraints dH) and ([5]). We do 
this by choosing a sequence of control inputs u* for each participant i that can vary with 5. 
We wish to choose the best causal response to prediction errors, a policy u* = vrj((5), where 
to be chosen before the error is known. "Causal" means that can 
depend only on the measurements of 5o, 5i, • • • , Sm, 6m+i- That is, we assume that 5m+i, the 
sub- vector of 5 pertaining to time m + 1, is revealed just before input u!^ is applied. Obviously, 
a dependence on any of 5m+2, ■ ■ ■ , ^T-i would violate causality because would be a function 
of information not yet available at time m. 

Substituting u* = 7rj(5) into the state update equation ([1]) and eliminating x\ we obtain 



the following finite horizon optimization problem: 



min y^^E[Ji{Aix'Q + BiTTi{6),7ri{6))] 

Causal vr,- ^— ' 
i=l 



s.t. ^n + Gid + Ci{Aixi + BiTTi{6)) =0,V5e A, 



i=l 



MS) 

6 



eZi,\/5 e A. 



(6a) 
(6b) 
(6c) 

(6d) 



Constraints (|6bp and ([6c]l are the expanded forms of (jH) and ([5]) after substituting definitions 
of X* and u* , so that the optimization is expressed only in terms of TTj . 

This problem is intractable due to the wide variety of candidate functions tTj that could 
satisfy the constraints. We therefore restrict ourselves from now on exclusively to policies of 
the affine form 

= + ei , (7) 



T-lJ 



SO that participant i's power schedule u' is parameterized by a nominal schedule = [bq . . . & 
plus a linear variation Di with future prediction errors. In order for the use of future distur- 
bances to be causal, Di takes the lower-triangular form 



[A]o,o 

[A]i,o [A]i,i 







\D 



jjT-1,0 



\D 



jjT-l,T-2 [AjT-l,T-l 



where [A]/,m G M^^^* is the response of input ul to error 6m+i- The presence of non-zero 
elements above the diagonal would violate causality, because the values of as-yet-unknown 
errors would contribute to the control rule. 

Constraint ()6dp and the causality requirement are rewritten for compactness as a set of 
admissible policies {Di,e.i) parameterized by the current state Xq: 



(A,ei 



[A]/,™ = 0, 

Axl + Bi{Di5 + e^ 
Di6 + Ci 



Mm > / 



GZi,V5 e A 



This leads to the fohowing rewriting of problem ^ in terms of nominal schedules {ei}^^-^ and 
policy matrices {Di}f^^: 

min Ji {xq ,Di,ei) (8a) 

(A,eOeJ-i(xj,) ^ 



s.t. J] rj + Gi6 + Ci(Ai4 + Bi{Di6 + e^)) = , 



1=1 



V(5 G A , (8b) 



Tiin + G^5 + C,(^ix^ + Bi{D,6 + e^))) < p , 



i=l 



V<5 G A . (8c) 



The objective function has been redefined as Ji{xQ, Di,ei) := E[Jj(x*,u*)] to reflect its 
dependence on Xq, Di, and ej. 

The modelling assumption of a positive semidefinite quadratic form ([2]) for Jj(x*, u*) allows 
the expectation over 5 to be minimized straightforwardly in comparison to the arbitrary case, 
since only the moments E[(5] and E[55'] are needed. Substitution from ([2]) gives the following 
representation of the objective, which is convex: 

Ji{xl, Di, a) = E[Ji(^ix'o + Bi{Di5 + ei),Di6 + e^)] 

= ff'{Aix\, + B,e^) + ffa + xyA[HfB,e^ 
+ \xl'A',HfA,xl + yMHfB, + Hf)e, 
+ {frB^ + fr + xl/A^H^Bi 

+ e[B[HfB, + e[Hf)Dim 
+ \{D[{B[HfBi + Fr)A,E[M']) + a 

Since the expected prediction error E[5] can always be set to zero (the reference predictions 
Tj would always be chosen with E[(5] = 0), the second-to-last term in general cancels. Note 
that although the constant term makes no difference to the solutions Di and Cj, it is needed 
later in order to compare the costs of different approaches. 

3.2 Equivalent tractable reformulation 

Problem ([8]) cannot be solved directly because constraints (|8bp and ([8c]) . as well as the defini- 
tion of J-^j(xQ), apply for all 5 G A and are therefore the intersection of an infinite number of 
constraints. To obtain a numerical solution they must be written in an equivalent finite form. 

Since A is assumed to have an interior containing the origin, i.e. no redundant components 
exist in 5, it can be shown that constraint ()8bp is satisfied if and only if the following conditions 



hold: 



^(ri + CiAixl + CiBiCi) = , 

Np 

J2iGi + QBiDi) = 0. 



(9a) 
(9b) 



Constraint ([8c|) and the sets J^i^Xq) can be written using a result due to Guslitser and 
others |14p i5j. and employed in [19] to solve a similar robust control problem to this one. 

Recalling that A = {5 | S'J < h}, the following equivalences hold for ([8c]l . An extra 
matrix variable Z is introduced for the last equivalence, which uses strong duality in linear 
programming. 

Np 



Y,^i{ri + G,5 + C,{A,xl + B,{D,5 + e,))) < V5 G A 

t 



N„ 



Nr, 



1=1 



<p-^ri{ri + CiAixi) 



i=l 



t 



BZ-.Z'h+Zfli T.CiBie, < p -Zfli ^iiri + QAixi), 
Ei=i ^iiGi + CiBiDi) = Z'S, and Z > element -wise. 

Similarly, it can be shown that the sets ^-^(xq) can be rewritten in finite form as follows, 
starting from definition ^ and introducing extra matrix variables Yi of appropriate dimension: 



{Di,ei 



[A]z,m = 0, Vm > Z 

3y, > 0: {T,B^ + Ui)D^ + Vi = ¥(8, 

TiAix'o + {TiB, + Ui)ei + y//i < Wi 



These changes lead to the following tractable representation of optimization ([8]): 




(10a) 



s.t. ^(ri + CiAixl + CiBiei) = , 



(10b) 




(10c) 



1=1 




(lOd) 



i=l i=l 



Y,^iiGi + CiBiDi) = Z'S. 



(lOe) 



i=l 



In summary, (jlObp states that the nominal schedules of power output changes Cj should 
track the base prediction; (jlOcp states that the rules Di used by the participants should 
together track any error vector 5 € A; (jlOdp and (jlOep ensure that line current limits should 
not be exceeded for any 5 G A. 

After solving (|10p the state and input trajectories x* and u* for a particular prediction 
error 6 can be computed by substituting the solution {Di, Cj) back into ([7]) and ([1]). 

3.3 Computational issues 

Problem pop is a quadratic program, which in principle can be solved even where many 
thousands of variables and constraints are present. We now comment on the size and struc- 
tural properties of the problem. Each vector has T elements, each matrix Di has N^T'^ 
elements (neglecting the fact that some of these are constrained to be zero), and each ma- 
trix Yi (introduced by definition of the sets ^^^(xq)) has qli elements, recalling that q is the 
number of constraints defining A and li is the number of constraints defining Zi. The ma- 
trix Z has 2qLT elements. Therefore the total number of primal optimization variables is 
Np{T + NsT^ + ql) + 2qLT, where I = ^ ^t'l ^i- The main computational cost arises from 
matrices Di which together have NpN^T"^ elements. Clearly, the problem size already reaches 
the thousands for modest parameter choices. 

However for computational purposes the structure of the problem is as important as the 
size. The form of (|10p exhibits a convenient linear coupling between participants i = 1, . . . , Np 
in both the cost function and the constraint set, and in fact lends itself straightforwardly to 
solution via a large-scale distributed solution method, such as the recently-revived ADMM [28] . 





4 Properties of affine reserve policies 



We now highlight some fundamental properties of policy-based reserves. We use the word 
"policy" loosely to refer to the matrix Di (rather than the nominal component Cj) in what 
follows. 

Property 1. Assume that for all participants [x\_^-^]i = u\. Constraint (jlOcp is equivalent to 

e£iA = -e£iG.. 

Recall that from the model definition, 

CiB, = [It ^[10ixK_i)]] 

[Bi]i ••• " 

[AiB,]i [Bi]i ••. ; 

: ■■. ■■. 

_ [A^-^Bi]i ••• \A,Bi]i [Bi]i . 

The assumption, as described in Section [2711 is that the first state (the participant's net power 
infeed) is set directly by the input decision at the previous time step, and not as a function 
of current states. Therefore the first row of Ai contains only zeros, and [Bi\i = 1. From this 
it can easily be shown that CiBi = It, from which the property follows immediately. 

As an example, if a single source of uncertainty were present that maps directly to the 
power output of one participant, one would have YldZi = —It-, so that the diagonal com- 
ponents of Di sum to —1, and the off-diagonal entries cancel. 

Property 2. No freedom exists in the choice of policies unless more than one participant uses 
them. 

This property follows directly from Property [H and is to be expected. Consider the exam- 
ple of a single generator supplying an uncertain fluctuating load. The generator must track the 
power requirement at every time instant and there is no freedom to plan this task strategically. 
If one or more additional generators were to become available, tracking could be planned in 
many different ways according to the characteristics of the generators, thanks to the extra 
degrees of control freedom the additional generators would bring. 



Bi ••• 

AiBi Bi ■ ■ : 

: ■•. ••. 

Aj~^Bi ■ ■ ■ AiBi Bi 



Property 3. Assuming that the mappings Gi are memoryless, i.e. Gi = lT®Gi, the sum over 
all participants of a given sub-diagonal entry in the D-matrix must be zero, i.e. ~ 

0, yi>m. 

This is seemingly counter-intuitive, because it implies that all the "interesting" parts of the 
participants' policies should ultimately cancel each other out. However through the definition 
of the feasible set J^j(xQ), the choice of Di has a bearing on what nominal trajectories Cj are 
feasible for a given participant. Since both Di and Cj determine the real behaviour of each 
participant, operating costs may be affected significantly. 

5 Implementation of policy-based reserves 

Reserve provision in electrical networks is increasingly a priced, rather than regulated, com- 
modity [27]. In order for policy-based reserve provision to be compatible with liberalized 
electricity markets, there must be a clear indication to market participants of the value of 
offering these more elaborate services. That is to say, a price must be placed on the require- 
ment for participants to respond in certain ways to wind prediction errors. In this section, 
the Lagrangian of problem ()10p is used to show that efficient market prices exist for policies 
and that an exact analogue of the standard Locational Marginal Prices (LMPs), which arise 
in electricity markets from network congestion over a power system [23], can be defined for 
policy-based reserve products. We then discuss the important practical implications of trying 
to implement policy-based reserves. 

5.1 Efficient prices for policies 

A partial Lagrangian of problem (fTO|) (keeping the constraints {Di,ei) € J^iix^) and Z > 
but relaxing all others), L(Z, Di, ei, . . . , Dat^, eiv^. A, 11, i/, ^'), can be formed by introducing 
the multipliers A G for violation of constraint (fTOb]) . 11 G M^^^«^ for (fTOcl) . u £ for 



(fTndl) . and e M2'f^^><^'5 for (fTnel) . The Lagrangian thus has the form 
L(Z,L>i,ei,...,L>Arp,eArp,A,n,i/, = 

Jiix'o, Di, ei) + AM ^ ri + Qylixj, + CiS^. 



+ ^n,^G, + Q5,Ay (11) 

(Np Np 
Z'h + J2 ^^CiB^e^ -p + Y,^i{n + CiA^x},) 
i=l i=l 

+ l^, 5^ r,(G, + CiBiDi) - Z'S^ 

After some rearrangement, this partial Lagrangian can be rewritten in separable form as 
L{Z, Di,ei,. . . ,DNp,eNp,\Ii,v,'^) = 

Np 

J2iMxh,Di,ei) - A^e, - (H,, A)) (12) 
1=1 

+ u'Z'h - Z'S) + /(A, n, u, ^) 

where / is constant with respect to the primal variables Di, ei, and Z. The Lagrange mul- 
tipliers from optimization problems are commonly interpreted as prices, and here we have in 
effect defined two prices. The first is a nodal power price 

A, :=-B^Ci{X + nu), 

consisting of a global component depending on A and a local component (induced by any line 
congestion present) depending on z^. This agrees with the standard derivation of locational 
marginal prices (LMPs) in optimal power fiow theory. The second is a matrix of reserve policy 
prices 

Ui :=-i?,'cKn + r:^), 

whose entries are the marginal value of each element of D^. Note that the elements of Ilj 
above the main diagonal are not required, since these are related to entries that would anyway 
result in non-causal responses to prediction errors, were they to be different from zero. 

As discussed in Section [H BiCi = / in most cases, so that usually Aj = — (A + T'-v) 
and Ilj = — (n + r^^*). The minus signs are a result of the sign convention used to write 
constraints (jlObp and (jlOcp . The identical form of Aj and Hj suggests that optimal reserve 
prices for policies exhibit the same locational variation as LMPs. 

Because the Lagrangian is separable, the terms X'^Ci and {Ili,Di) can be identified as 
the efficient payments to each participant i for committing to nominal plan and reserve 



policy Di. The payments represent the money transfers that would result from a market 
mechanism that solved problem (jlOp efficiently. Under such a scheme, the expected profit 
made by participant i would be equal to — Jj(a;Q, Dj, Cj) + A-Cj + (Ilj, i^j). 



5.2 Policy-based reserve products 

A natural question is what physical commitment a participant makes when implementing a 
given policy (-Dj,ej). Clearly the nominal part Cj is the schedule participant i will follow if 
predictions turn out to have been made with perfect accuracy, i.e. (5 = 0. The interpretation 
of Di is more subtle and can be termed in two ways, which will be described here. In this 
subsection the analysis is given for a single uncertainty source, i.e. Ns = 1; the case for Ns > 1 
is analogous. 

Firstly, each row I of Di can be read as the rule the participant must follow to construct 
its input which for the realized errors 5q^...,6i determines the power it injects at time 



m,=0 

Secondly, each column m Di (reading down from the diagonal) is, in control terminology, 
the planned impulse response g{5m+i) of participant i to a unit prediction error m steps from 
the current time (recalling that 5m+i is revealed just before u]^ is applied): 



Example product 

Consider the example of the decision rule governing which sets the power participant i 
supplies during step / + 1 of the planning horizon. For given choices of [-Di];,05 to 
[-Djji^i agreed, the payment (using the price notation from Section [5.ip made to participant i 
for the reserve service would be X]m=oP*km[-Cj]«,m and the payment for scheduled power [cj]; 
would be [Aj]j[ej];. An analogous product could be sold based on the column-wise reading of 
Di, in which case the participant would be selling an a posteriori response to errors. 

Existing reserve mechanisms could be modelled by matrices Di for which only the main 
diagonal is populated. In secondary reserves provided by conventional generators, feedback 
controllers use the frequency deviation, in the form of the Area Control Error (ACE) signal, to 
adjust their power outputs to follow load mismatches. For a unit deviation from the net load 
reference, each generator will end up with an offset that depends on its controller parameters. 
For such a unit load deviation at time /, this offset is in our notation exactly the matrix entry 





Therefore, a simple way of comparing policy-based reserves to existing mechanisms is to 
restrict the structure of the participants' D-matrices accordingly and inspect the results. This 
is demonstrated in Section O 

5.3 Real time operation 

For continuous operation of the power system, it is necessary to choose new policies period- 
ically, because the policies only apply for the following T steps at the time when they are 
chosen. Furthermore the optimal policies are a function of the current state of the system, 
and it may be attractive to choose new policies early in the light of new forecast information. 
Therefore a systematic way of choosing policies repeatedly is required. Two possible schemes 
for such "closed-loop" operation are: 

• Batchwise: Set policies for a horizon of T steps, let all T steps play out, then choose a 
new batch of policies once the current ones have expired. 

• Receding horizon: Let only the first step play out, then immediately update the 
policies for the next T steps. 

Clearly, a middle road between these two options also exists, namely letting some of the T 
steps play out, then discarding the remainder and choosing new policies. 

The receding horizon scheme presents another choice — whether to honour previous poli- 
cies in the new choices for Di (respectively Cj). This would be done by shifting the matrices 
(resp. vectors) upward and leftward (resp. upward) by one element, then optimizing over the 
bottom row of elements (resp. element) to obtain new policies. These new policies would 
be feasible with respect to J^i(xg) as well as the global constraints ()10bp - ()10ep under mild 
assumptions on A (for conciseness this is stated without proof here). 

An alternative to this is to reject the previous policies and choose fresh values for Di and 
Cj at every step. This option is attractive because it allows new (presumably better) policies 
to be chosen in the light of new information as soon as it becomes available. However its 
drawback is that due to the repeated renewal process, only the (0, 0) block of the D-matrix 
ever gets used, and no policies beyond the first row ever appear to be implemented. This 
conflicts with the ideas developed in Section 15.21 in that contracts for policy-based reserves 
would be agreed but then never called on. It is important to note, however, that the optimal 
values of [e^Jo and [I?i]o,o would be affected strongly by the costs and constraints modelled for 
steps 1 to T — 1, even though those later steps would never be realized. This means that under 
this scheme, any cost savings from employing reserve policies are ultimately to be found in 
the more intelligent choice of [ej]o and [Z)j]o,o- This effect is reported in Section [6l 

Rejecting previous policies at each time step makes choosing correct payments to market 
participants difficult. One way of overcoming this would be instead to define the new policies 



as deviations from shifted versions of the existing previous pohcies, with payments settled 
additively. Under such a scheme there would be T payments for the reserve action undertaken 
at a given time, resulting from the superposition of policy choices made over the last T steps. 

6 Numerical case study 

We now apply policy-based reserves to a standard test network that has been adapted to 
include a large share of wind power. Policies are recomputed in a receding horizon fashion 
(the second method described in Section 15. 3p over a three-day simulated period in order to 
assess the reserve costs incurred. The test is repeated 100 times with different realizations of 
the random wind infeed. The effects of restrictions on the structure of the matrices Di on 
the cost of reserves are reported, leading to observations on the cost savings facilitated by 
recourse. The optimization problems have been solved directly using CPLEX. 

The network used is a modification of the 39 bus network described in Appendix A of 
[29], and shown in Fig. [TJ This network contains 7 thermal generators, 2 storage units, 19 
loads, and 3 wind farms (which replace 3 of the original 10 generators). We assume that the 
generators represent the plants that have been selected for use via an earlier unit commitment 
decision. The parameters in terms of the definitions in Section 12.11 are described in Table 
[TJ The daily pattern of load variation shown with the thick black line in Fig. [3] was taken 
from data for total UK national electrical consumption on 14th September 2012 (available 
at http://www.bmreports.com), normalized to the size of each load modelled. Peak load is 
around 6 GW, and wind power provides 3 GW at maximum. Line flows from bus 16 to 15 and 
from bus 16 to 17 were restricted to 1000 MW. The load sizes Pnom in Table [1] are the nominal 
values described in [29]. The simulation period was three days, divided into 288 fifteen-minute 
time steps. The horizon length was T = 8 steps. 

6.1 Uncertainty model 

Uncertainties in the system were assumed to originate only in the random wind power avail- 
ability (loads were assumed to be predicted exactly). The wind farm output was driven by 
the following first order random process model with saturation: 

Qk+i = min{max{gmin, Qk + A}, gmax} , 

where qk G denotes the state of the uncertainty model, and (3k ~ -^(0, S) is sampled at 
each step k from a multivariate normal distribution with variance S € W^^^^^, truncated to 
bounds A^/3fc < 6^. 

Defining q := [q'l ■ ■ ■ q'j']' as the random future evolution of q from current state qo, the 
nominal predictions of wind farm power output rj were mapped linearly from E[q], so that 




Figure 1: 39 bus test network from [29], with wind infeed replacing thermal generators at 
nodes 32, 33, and 34, and with added storage units at nodes 1 and 28. 

rj = GjE[g], where Gi is the same matrix as that described in Section l2.ll The prediction 
error was defined as 6 := g — E[(7] so that E[5] =0. The prediction error covariance is then 
E[(5(5'] = E[((7 — E[g])((7 — E[(7])']. This was supplied together with the references r, and the 
current system state Xq as inputs to optimization problem (jlOp . 

The uncertainty set A was recomputed as a function of the current state, to reflect the 
fact that prediction errors are bounded relative to the nominal predictions by the wind farm 
power output limits and by the bounds assumed on the change in wind power availability 
at each step. Its new parameters 5 and h were then supplied to (fTO|) . In addition, at every 
simulation time step, an aggregation of 20,000 Monte Carlo runs was used to produce T-step 
nominal predictions rj for each wind farm, as well as estimates of E[5(5'], the variance of the 
error relative to the nominal value. 

The parameters defining the uncertainty model in this case study are given in Table [2j The 
three wind farms are driven by two temporally and spatially correlated sources of uncertainty. 
The matrices Gi for the wind farms are of the form Gi = lT®Gi, and the matrices Gi are given 
in Table [T] The output of the wind farm at node 34 is the mean of the two at nodes 32 and 
33. The purpose of this formulation is to show that random power flows can be driven by a 
process with dimension lower than the number of nodes affected, either to save computational 
effort, or because the uncertainty is anyway difficult to model. 

An example time series qk generated in this way is shown in Fig. [21 

6.2 Comparing reserve costs 

For each of the 100 wind realization tests, three parallel models of the power system were 
driven by the same realizations of the random wind model described above. The observed 
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Figure 2: Example output of correlated uncertainty model g^, which drives the wind farm 
outputs in the case study. The solid line represents [qk]i and the dotted line [qk]2, for the 
-j^QQth realization used. 

operating costs under receding horizon control over the simulation period were then compared 
under the three schemes, which were defined as follows: 

1. Prescient case: Disturbances are known at the time the finite-horizon optimization is 
carried out. Matrices Di are therefore not needed, and nominal schedules track the 
power reference perfectly. This model, which results in the best attainable receding- 
horizon cost, is used as a point of comparison for the other two. 

2. Flexible-rate reserves: [-Dj]z,m = for / ^ m, for all i. This represents the best possible 
response to uncertainty without using recourse, and the optimization is over the elements 
of ei and the diagonal parts of Di.. 

3. Policy-based reserves: [Dili^^ = for / < m, for all i. This allows full use of the extra 
information that will be available at each time step when the reserve is deployed. 

The total operation cost was measured from the state and input values x^{t) and n*(t) 
realized by the elastic participants at each time step t over the simulation period. 



The inputs are indexed by {t — 1) because state x(0) is given whereas input n(0) must be 
chosen and determines and so on. 

An example of the power output traces for the generators is given in Fig. [3l The cost results 
are shown in Table O A cost of reserves is defined for Schemes 2 and 3 as the operation cost 
experienced minus the prescient cost (from Scheme 1). This represents the cost incurred in 
order to accommodate the uncertain wind infeed. Across the 100 runs, the cost of reserves 
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decreased by an average of 34.5 % for full policies (Scheme 3) with respect to the best possible 
non-recourse reserve scheduler (Scheme 2). Average costs under Scheme 1 were 4.318 x lO''' 
for the three-day test, and reserves were found to incur up to 1 % of additional total power 
generation costs in the tests. 

Another test was carried out allowing only the diagonal and immediate sub diagonal entries 
in the D-matrices to differ from zero, representing a "mixture" of Schemes 2 and 3. The 
average savings in the cost of reserves, measured in the same way as for Scheme 3, were 32.4 
%, indicating that most of the 34.5 % savings can be obtained from using just a two-step 
policy rather than a full policy (in this case 8 steps). Whether this is always true, however, is 
likely to depend strongly on the parameters of the participants connected to the grid. 



Trace of generator and storage output 




Time step 

Figure 3: Power output traces for the 100**^ wind realization test under Scheme 3. Generator 
outputs are plotted in stacked form, and the total power output of the two storage units is 
plotted using the thin black line relative to the top of the stack. The wind power injection is 
the difference between the total load (bold black line) and the sum of storage and generator 
infeeds (thin black line). 

7 Conclusions 

Motivated by the need to incorporate intermittent renewable energy sources into power sys- 
tems as intelligently as possible, this paper introduced the idea of policy-based reserves. The 
approach relies on solving a robust optimization problem to find efficient responses to errors 
in the prediction of uncertain load or supply, over a finite horizon. Interpretations were de- 
veloped for policy-based reserves as products that could be bought or sold on power markets, 
with associated prices. A case study demonstrated the operating principles and benefits of 
the approach, and insights were gained into the important properties of policy-based reserves. 

The cost of reserves measured in this paper were determined in large part by the correlation 
of the prediction errors. This is evident from the term ^{D'-{B'-HfBi + Hy')Di,E[d6']) in the 
cost function. Correlation in space and time results in non-zero off-diagonal entries in E[5(5'], 



which will change the contributions of individual entries of Di to the expected costs. It is not 
surprising, then, given the strong correlation of prediction errors found in the numerical case 
study, that for a given finite-horizon optimization the costs were reduced when recourse was 
allowed. 

In receding horizon operation, new power schedules are found at every time step. This 
can be viewed as a form of recourse that is not taken into account for each finite horizon 
optimization. This leads to the intuition, borne out by the numerical results, that some of 
the apparent value of optimizing over full lower-triangular policies (rather than diagonal ones) 
should disappear once the optimization repeated in receding horizon fashion. However, as 
reported, the cost reduction can still be significant. It was also observed that most of the cost 
reduction can be gained using only one or two subdiagonal rows in the D-matrices. 

Lastly, the interpretation of the reserve cost defined in the case study is important. Scheme 
1 represents a lower bound on the cost achievable under any possible receding horizon scheme, 
however it may not be possible to approach this value using any causal policy. It is unknown 
whether the afhne policies used are close to being the best possible reserve mechanism for this 
system model, or whether some other policies 7rj(5) could be found that would result in even 
lower reserve costs. Recent related results on finding performance bounds for affine control 
policies [30] may provide some insight into this issue. 
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Table 1: Parameters of elastic participants 
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Storage units, i = 8,9: 
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Wind farms [Node: Gi 



32: [2 0], 33: [1 1], 3^: [0 2] 



Loads [Node: pnom] 3: 322, 4: 500, 7: 233.8, 8: 522, 12: 7.5, 
15: 320, 16: 329, 18: 158, 20: 628, 21: 274, 23: 247.5, 24 : 308.6, 
25: 224, 26: 139, 27: 281, 28: 206, 283.5, 31: 9.2, 5f; 1104 
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e 3: Cost comparison for different structural restrictions on 
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